In this short course we will introduce some basic concepts of data oriented programming, both in R and Python.
A short list of the topics that will get covered is:
- Methods of programming
- Basic of programming
- Data structures and databases
- Bootstrap techniques
This concept will be discussed separaterly in the two lessons, one using R and the other using python.
All the concepts are applicable in both languages. One can use one language or the other, or a mixture or both.
Methods of programming¶
- shell vs scripting
- dividing a program in pipelines
- data formats, text vs binary
- complete automation and configuration
- random data generation, random seeds and reproducibility
- readibility and longevity of a program
Basic of programming¶
- cycles and iterations
- map, filter, reduce
- vectorization
- functions
- DRY principle (Don't Repeat Yourself)
- functions and routines
- isolation principle and functional programming
Data structures and databases¶
- databases and dataframes
- tidy data and normalization principles
- pivot, join, stacking and SQL
- split-apply-combine pattern (groupby)
Bootstrap¶
- visualization with tidy data
- uncertainty in report and visualization
- uncertainty for model checking
- data generation, bootstrap, subsampling
programming and coding¶
This two concept are often mixed up. We could use this two informal descriptions:
Programming is about thinking how to formalize a certain sequence of operation.
Coding is writing a computer program that actually perfom the operations we have in mind.
While coding is strictly related to a specific programming language, but most of the programming skills are identical among different languages.
As a matter of fact, one does not even need a computer to be able to "program"
how do you prepare coffee?¶
Try describing the steps necessary to prepare a coffee.
pointers¶
- are you describing a moka or a capsule machines?
- how expert is the other person?
- is the coffee already available?
shell Vs scripting¶
Even if the language is the same, the goal is completely different between the two.
Shell programming is useful for explotatory analysis, trying to understand what it's going on in your data
Script programming is useful to perform your analysis in a formal and reproducible way.
reproducibility is key
- six months from now, you have to redo a figure for a paper because you realize you put the wrong labels
- four years from now, they will ask you how you did that analysis, in details
- tomorrow you decide to change a passage in between
pipelines¶
divide your program in meaningful, separated scripts.
- gather raw data
- clean the raw data into tidy data
- transform the data for analysis
- do the analysis
- visualize them
each step should have its own script, that take inputs and gives output used by the following program
you need to repeat only the few last steps, not the whole thing, and can do that with the press of a button
data storage¶
Never, never, never, delete or modify the raw data!
When possible, choose text over binary format. Computer space is cheap, understanding the data is not
Metadata is the new gold
Make sure your data is always computer readable!
metadata¶
- When was the data collected?
- by whom?
- how?
- there were any problems?
- what each variable means
- what are reasonable values for it
In the name of everything that is holy, never put the metadata as inplicit information. This includes:
- color of an excel cell
- test formatting in excel
- excel comments
- case of the text
- arcane special values (-9 for missing, -99 for unvailable, etc...)
If it is necessary to update the values with new informations, keep an update table on the side and incorporate it into the data cleaning process.
| table name | key value | variable name | note |
|---|---|---|---|
| phenotypical | John Black | weight | data is missing because the patient refused |
| biochemical | Jane Doe | cholesterol level | the vial was dropped and the sample lost |
| medical_history | Dave White | age of menarche | data is not applicable |
Your code to load and analyze the data is part of the metadata!
A report is nice, but it might be inaccurate.
The code describe exactly what has been done to the data.
Your environment configuration (OS version, R version, packages version) is part of the metadata as well.
Different versions of the same libraries could give different results!
sessionInfo()
packageVersion("tidyverse")
Today we will discuss the basic of data exploratory analysis with R.
Wednesday we will discuss more data wrangling and model checking using Python.
Today I will show the similarities between python and R syntax for data exploration, and try to show some differencies.
programming concepts¶
There are really only 4 basic concepts to programming; anything else, in all languages, is basically just rehashing of these ideas.
- Conditional execution
- functions and routines
- iteration
- data structure(s)
conditional execution¶
do something only if a certain condition is met
x <- -5
if(x > 0) print("Non-negative number") else print("Negative number")
The if and else clauses can contain several statements, just be careful to keep the else on the same line as the closing brace
x <- 5
if(x > 0){
print("Positive number")
} else {
print("negative number")
}
x <- 5
y <- ifelse(x>0, 1, -1)
y
Python has a very similar look, but with less clutter:
x = 5
if x>0:
print("positive number")
else:
print("negative number")
y = 1 if x>0 else -1
functions and routines¶
From the language perspective of both python and R they are the same thing, but from a concept perspective, they are very different.
- routines do something to the system, and usually don't return anything interesting (like printing).
- functions create new objects from the inputs, but don't modify the input themselves
In both languages they are referred as functions, and can do both, but try to know which is which and plan accordingly
function_sign <- function(x_value){
sign <- ifelse(x_value>0, 1, -1)
return(sign)
}
function_sign(3)
in python, again, they are very similar
def function_sign(x_value):
sign = 1 if x_value>0 else -1
return sign
function_sign(-3)
R functions give some advanced functionalities, like partial recognition of the arguments.
Please, avoid relying on them, only evil can come out of it...
function_sign(x_value=-3) == function_sign(x=-3)
An inportant concept related to functions is the concept of scoping.
When we execute a function we can use variables that are reserved to the function.
This allow to write simpler code.
In R the scoping is the major difference between assigning with = and <-.
the easy solution is to always use <-
my_adder <- function(x){
g <- 3
return(x + y + g)
}
y = 4
print(my_adder(1))
print(g)
DRY principle¶
DRY: Do not Repeat Yourself.
Repeating code (or copy pasting it) can lead to sever errors:
- dependencies on surrounding code
- forgetting to change a variable name
- visual clutter
Functions are the primary way of implementing DRY principles.
The rule of thumb is that if you are going to do it more than twice, write a function instead... Sometime even if you are going to do it once!
x <- rexp(100, rate = 1)
x_med <- median(x)
x_min = min(x)
x_max = max(x)
x_top = x_max - x_med
x_bot = x_med - x_min
x_ratio = x_top/x_bot
x_ratio
This code is obscure and pollute the namespace!
interquantile_asymmetry <- function(data){
x_med <- median(data)
x_min = min(data)
x_max = max(data)
x_top = x_max - x_med
x_bot = x_med - x_min
x_ratio = x_top/x_bot
return(x_ratio)
}
interquantile_asymmetry(x)
iterations¶
Repeating an operation over and over, tipically with a different input each time, or with an exit condition.
Typically they underpin the usage of words like: while, for/for each, until.
from datacamp
years = c(2010, 2011, 2012)
for (year in years){
print(paste("The year is", year))
}
in python:
years = [2010, 2011, 2012]
for year in years:
print("The year is", year)
x <- 0
while (x<3){
print(paste("step", x))
x <- x+1
}
in python
x = 0
while x<3:
print("step", x)
x = x+1
data structures¶
This are containers that holds the informations, and make more or less easy to retrieve and manipulate this information.
Programming can be described as creating and modifying data structures with the previous concept until the desired result is obtained.
There are a great number of data structures, and more get introduced in a language by libraries.
The most common and important one for our goals are tabular data.
tabular data¶
Tabular data can be imagined as a single sheet of an Excel Spreadsheet, or a database table.
There are columns which describe different kinds of informations, and rows that contains the various informations for each observation (tuples)
in R they are implemented by the native data.frame or the more advanced tibble from dplyr.
In Python they are implemented as DataFrames from the pandas library
vectorization¶
vectorization is the ability of a programming language to repeat a specific operation on all the elements of a data structure.
It is a form of inplicit iteration, with the possibility of inplicit conditionals.
It is an extremely important capability for data analysis, as it allow to express complex idea in an easy way
x <- c(1, 2, 3, 4)
x + 1
for (xi in x){
print(xi + 1)
}
x[x < 3] * 2
for (xi in x){
if (xi < 3) print(xi * 2)
}
csv and tsv files¶
the most common way one can obtain tabular data is by obtaining a csv or a tsv (comma separated values and tabular separated values), a text file encoding the table data explicitly.
My personal preference goes toward tsv, as they are more human readable:
a b c
170 M 1983
180 F 1972
compared to the equivalent csv:
a,b,c
170,M,1983
180,F,1972For this course you should have been given few csv files:
- UPLC_Plasma_Clinic.txt
- UPLC_Plasma_ExpDes.txt
- UPLC_Plasma_QC.txt
- UPLC_Plasma_Raw.txt
notice that, even if the file name does not end with csv, they are perfectly valid anyway.
A csv is just a normal text file, formatted appropriatedly!
When in doubt, open it and check the formatting...
import¶
We have functions that read the data files and convert them directly into data.frames.
Never try do to this by hand, unless you are very certain of what you are doing, there are just too many pitfall
This is the first step to get your data in the right format.
This functions have several nice abilities, including automatic data type recognition: dates get recognized as dates, numbers as numbers.
When they fail, we need programming abilities to supplement their work.
UPLC_Plasma_Clinic <- read.csv("data/UPLC_Plasma_Clinic.txt",
header=TRUE, sep='\t', dec=',')
head(UPLC_Plasma_Clinic)
the equivalent code in python would be:
import pandas as pd
UPLC_Plasma_Clinic = pd.read_csv("data/UPLC_Plasma_Clinic.txt",
sep='\t',
dec=',',
index_col='Sample')
UPLC_Plasma_Clinic.head()
tidy¶
Format as tidy data.
Tidy data are a sane way of representing you data in tabular form (correspond to the third normal form of database tables)
Each columns is a variable you want to observe, each row is an observation unit.
for R, the reference liset of libraries to work with tidy data is called tidyverse
library("tidyverse")
There is a specific function for loading data using the tidyverse function read_csv (or read_tsv or read_delim).
They will load it into a tibble, a data structure similar to a basic R data.frame, but with some adjustment.
tibbles are a more modern structure, and when possible I suggest using them.
UPLC_Plasma_Clinic <- read_tsv("data/UPLC_Plasma_Clinic.txt", locale=locale(decimal_mark = ","))
head(UPLC_Plasma_Clinic)
There are various operations, defined as single verbs functions, that take the data.frame as input and send out another data.frame as output.
This operations can be linked as pipes (even if I personally would recommend against it)
- Pick observations by their values (
filter). - Reorder the rows (
arrange). - Pick variables by their names (
select). - Create new variables with functions of existing variables (
mutate). - Collapse many values down to a single summary (
summarise). - and so on...
The first argument is a data.frame.
The subsequent arguments describe what to do with the data frame, using the variable names (without quotes).
The result is a new data frame, that can be passed to other functions in a similar fashion
partial_data <- filter(UPLC_Plasma_Clinic, CaseControl == 1, Sex=='male')
head(partial_data)
in python there would be expressed as
table[table.column_name == some_value]
table[table['column_name'] == some_value]
or
table.query('column_name == some_value')
All these functions can be used in a functional way:
query = pd.DataFrame.query
query(table, 'column_name == some_value')
In general R manipulation functions looks more consistent than Python's one, as pandas's DataFrame try to keep a behavior consistent (where possible) with standard python.
beware of floating point numbers when doing comparisons!
sqrt(2) ^ 2 == 2
1/49 * 49 == 1
this is a problem of how computers represents numbers, not of R.
One would obtain the same results in R, Python, Matlab, C or most programming languages.
Could be solved with symbolic algebra, but it's a whole different topic (both R and Python support symbolic algebra witht external libraries).
Make it tidy (check for problems)¶
In our case the data that we have is already in a tidy format.
at least apparently.
We should always check to make sure that the assumptions that we hold are correct.
This assumptions should be encoded in the script, so that they are maintained over time and these assumptions can ben checked again explicitly
nrow(UPLC_Plasma_Clinic)
length(UPLC_Plasma_Clinic)
stopifnot(length(UPLC_Plasma_Clinic)==4)
stopifnot(nrow(UPLC_Plasma_Clinic)==1235)
unique(UPLC_Plasma_Clinic$Sex)
stopifnot(length(unique(UPLC_Plasma_Clinic$Sex))==2)
stopifnot(min(UPLC_Plasma_Clinic$Age)>0)
unique(UPLC_Plasma_Clinic$CaseControl)
stopifnot(length(unique(UPLC_Plasma_Clinic$CaseControl))==2)
display¶
Visualizing data is extremely important because it can create a more intuitive understanding of what the data looks like.
the most common library used for visualization in R is called ggplot, an implementation of the concept described in the grammar of graphic review.
ggplot make it very easy to create and manipulate graphs in different ways
ggplot(data=UPLC_Plasma_Clinic, aes(Age)) + geom_histogram()
ggplot(data=UPLC_Plasma_Clinic, aes(Age)) + geom_histogram() + facet_grid(Sex ~ CaseControl)
in python si potrebbe fare lo stesso con una notazione quasi identica con la libreria plotnine:
from plotnine import ggplot, geom_histogram, facet_grid, aes
ggplot(UPLC_Plasma_Clinic, aes('Age')) + geom_histogram() + facet_grid('Sex ~ CaseControl')
oppure usando seaborn, ma con una sintassi differente:
import seaborn as sbn
fg = sbn.FacetGrid(data=UPLC_Plasma_Clinic, row='CaseControl', col='Sex')
fg.map(plt.hist, "Age")
python is often a little more cluncky for the display, requiring more time to obtain a good looking one, but gives more power and control to the programmer to modify the graph as they seems more fit
transform¶
now we load the data from the experimental design, detailing in which plate each experiment has been done, and want to check if there was an unbalance of gender or age in each plate.
UPLC_Plasma_ExpDes <- read.csv("data/UPLC_Plasma_ExpDes.txt",
header=TRUE, sep='\t', dec=',')
head(UPLC_Plasma_ExpDes)
nrow(UPLC_Plasma_ExpDes)
nrow(UPLC_Plasma_Clinic)
wait...
let's re-read the details of the data:
Replicated standards have "stand" label and duplicated samples have "_D" label
so, if we ignore the replicated, we should have the same index, right?
library(stringr)
partial_data = (UPLC_Plasma_ExpDes
%>% filter(!str_detect(Sample, 'stand'))
%>% filter(!str_detect(Sample, '_D')))
head(partial_data)
print(c(nrow(partial_data), nrow(UPLC_Plasma_Clinic)))
print(length(unique(UPLC_Plasma_Clinic$Sample)))
print(length(unique(partial_data$Sample)))
We need to start working with out transform to see what is going on.
There are two main families of transformations that one needs to understand:
- grouping
- joining
- pivoting
grouping divide the dataset in sub dataset based on certain properties, apply an operation and merge the result together (such as calculating the average age for each Sex).
joining merge two different tables.
pivoting is used to transform something like a tidy table in a more squared table. The inverse operation is called melt.
The merging can be applied in several diffent ways, but the great four are:
- inner join
- outer join
- left join
- right join
from R for data science
there are other transformation that one can apply, such as creating new columns from old ones, removing or changing the values of various columns, but these don't require a detailed analysis.
anti_join(partial_data, UPLC_Plasma_Clinic, by='Sample')
joined = left_join(UPLC_Plasma_Clinic, UPLC_Plasma_ExpDes, by='Sample')
head(joined)
grouping¶
grouped_by_gender_and_plate <- group_by(joined, Sex, Plate)
summarized = summarise(grouped_by_gender_and_plate,
Age_mean = mean(Age, na.rm = TRUE),
)
summarized
pivoting¶
gender_age_per_cell = spread(summarized, key = Sex, value = Age_mean)
gender_age_per_cell
gender_age_per_cell['female'] - gender_age_per_cell['male']
to find out if it is significant, we have to verify if this difference is big or not compared to the uncertainties of the averages.
try to make a plot visualizing this!
hints¶
- potete fare un scatterplot con gli errori con le funzioni
geom_pointegeom_errorbar - per fare il grafico vi serve un singolo dataframe in versione tidy con tutte le informazioni
grouped_by_gender_and_plate <- group_by(joined, Sex, Plate)
gender_age_per_cell = summarise(grouped_by_gender_and_plate,
Age_mean = mean(Age, na.rm = TRUE))
grouped_by_gender_and_plate <- group_by(joined, Sex, Plate)
gender_sem_per_cell = summarise(grouped_by_gender_and_plate,
Age_sem = sd(Age, na.rm = TRUE)/mean(Age, na.rm = TRUE))
joined_2 = inner_join(gender_age_per_cell, gender_sem_per_cell, by=c('Plate', 'Sex'), suffix=c('_mean', '_sem'))
head(joined_2)
(ggplot(joined_2, aes(x=Plate, y=Age_mean, color=Sex))
+ geom_point()
+ geom_errorbar(aes(ymin = Age_mean-2*Age_sem, ymax = Age_mean+2*Age_sem)))
grouped_by_gender_and_plate <- group_by(joined, Sex, Plate)
gender_age_per_cell = summarise(grouped_by_gender_and_plate,
Age_mean = mean(Age, na.rm = TRUE))
grouped_by_gender_and_plate <- group_by(joined, Sex, Plate)
gender_sem_per_cell = summarise(grouped_by_gender_and_plate,
Age_sem = sd(Age, na.rm = TRUE))
joined_3 = inner_join(gender_age_per_cell, gender_sem_per_cell, by=c('Plate', 'Sex'), suffix=c('_mean', '_sem'))
(ggplot(joined_3, aes(x=Plate, y=Age_mean, color=Sex))
+ geom_point()
+ geom_errorbar(aes(ymin = Age_mean-2*Age_sem, ymax = Age_mean+2*Age_sem)))
UPLC_Plasma_QC = read.csv("data/UPLC_Plasma_QC.txt", header=TRUE, sep='\t', dec=',')
head(UPLC_Plasma_QC)
head(anti_join(UPLC_Plasma_QC, UPLC_Plasma_Clinic, by='Sample'))